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ABSTRACT 

We have modeled the simultaneous first-order Fermi shock acceleration of protons, electrons, and 
helium nuclei by relativistic shocks. By parameterizing the particle diffusion, our steady-state Monte 
Carlo simulation allows us to follow particles from particle injection at non-relativistic thermal ener¬ 
gies to above PeV energies, including the nonlinear smoothing of the shock structure due to cosmic-ray 
(CR) backpressure. We observe the mass-to-charge {A/Z) enhancement effect believed to occur in 
efficient Fermi acceleration in non-relativistic shocks and we parameterize the transfer of ion energy to 
electrons seen in particle-in-cell (PIC) simulations. For a given set of environmental and model param¬ 
eters, the Monte Carlo simulation determines the absolute normalization of the particle distributions 
and the resulting synchrotron, inverse Compton, and pion-decay emission in a largely self-consistent 
manner. The simulation is flexible and can be readily used with a wide range of parameters typical 
of 7 -ray burst (GRB) afterglows. We describe some preliminary results for photon emission from 
shocks of different Lorentz factors and outline how the Monte Carlo simulation can be generalized 
and coupled to hydrodynamic simulations of GRB blast waves. We assume Bohrn diffusion for sim¬ 
plicity but emphasize that the nonlinear effects we describe stem mainly from an extended shock 
precursor where higher energy particles diffuse further upstream. Quantitative differences will occur 
with different diffusion models, particularly for the maximum CR energy and photon emission, but 
these nonlinear effects should be qualitatively similar as long as the scattering mean free path is an 
increasing function of momentum. 

Keywords: acceleration of particles — ISM: cosmic rays — gamma-ray bursts — magnetohydrody¬ 
namics (MHD) — shock waves — turbulence 


1. INTRODUCTION 

Efficient first-order Fermi shock acceleration (also 
called diffusive shock acceleration) is often suggested as a 
likely mechanism for converting the bulk kinetic energy 
of relativi stic plasma flows into individual particle en - 
ergy fe.g.. IBvkov fe Treumannll20TTl: IBvkov et al.l[20T2li . 
However, many aspects of particle acceleration in rela¬ 
tivistic shocks remain uncertain because of the inherent 
complexity of the process. The particle distributions are 
highly anisotropic and the magnetic turbulence, essential 
for acceleration to occur, must be self- generated and is 


extre m ely difficult to c h aracterize ( Lemoine_&_^elletiei 


2003 : iNiemiec fe Ostrowskil 120061: 


Lemoine fe Pelletiei 


2010 ; iRcvillc fc Belli 12014 iLemoine et al.ll2014lh These 
roadblocks can be overcome with particle-in-cell (PIC) 
simul ations and intensive work has been done in this are a 
(e.g., iNishikawa et al.l[20071 : iSironi fe Spitkovskvl 1201 111 . 
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However, current PIC simulations are computationally 
costly and have a limited dynamic range. The trans- 
relativistic regi me, which may be important for GRB 
afterglows fe.g.. lMeszarod[200ffi lAckermann et al.ll2013l) 
and some types of supernovae fe.g.. IChakraborti et al l 
12011 , is less well explored eit her analytically o r using 
PIC simulations (see, however. iCasse et al.ll2013lh 

In this paper we model the nonlinear acceleration of 
electrons and ions (protons and He 2+ ) at relativistic col¬ 
lisionless shocks using a Monte Carlo simul ation of first- 
order Fermi shock acceleration (e.g. , lEllison fc Doublel 
l2002t lEllison. Warren, fc Bvkovll2f) bill . The steady-state 
Monte Carlo simulation parameterizes magnetic turbu¬ 
lence generation and particle diffusion, important ap¬ 
proximations but ones that allow a large dynamic range; 
the simultaneous acceleration of ions and electrons (along 
with the radiation they produce); and a self-consistent 
determination of the shock structure. We believe this 
is the first attempt, apart from PIC simulations, to in¬ 
clude electrons self-consistently with ions in a nonlinear, 
relativistic, Fermi acceleration. 

While it is well known that nonlinear Fermi accel¬ 
eration should preferentially inject and accelerate high 
mass-to-charge particle s compared to protons in non- 
relativistic shocks (e.g ., lEllison et a1.lll98l1 lEichleHll984 
Hones fe Ellisonlll991[f . to our knowledge this process has 
not been investigated in relativistic or trans-relativistic 
shocks until now. The non-relativistic 11 A/Z" effect (A 


































2 


Warren, Ellison, Bykov, Lee 


is the mass in units of the proton mass m p and Z is the 
charge in units of the electron charge e) has been shown 
to be consistent with observati ons of diffuse ions a ccel- 
erated at the Earth bow shock (lEllison et alJll990lf and 
has been used to model the shock acceleration of inter¬ 
stellar gas and dust, matching important aspects of the 
galactic cosmic-rav abundances observed at Earth (e.g., 
Ellison et al.l 119971: IMever et al.l 119971 IMever fc Ellisonl 

19991: iRauch et aLfl201d: IB : nns ct al. loii 

The A/Z enhancement we model is a purely kine¬ 
matic effect. It depends on the following assumptions: (i) 
the acceleration process is efficient enough so the back¬ 
pressure from accelerated particles noticeably modifies 
the shock precursor, (ii) all particles have a scattering 
mean-free-path of the approximate form 

Amfp — ^mfp/mfp ? ( 1 ) 

where / m f p is an increasing function of local frame mo¬ 
mentum p, as is generally assumed, and (iii) the normal¬ 
ization parameter 77 m f p setting the “strength” of scatter¬ 
ing is similar for all particle species. 

The simplest assumption for particle diffusion is that 
/mfp equals the gyroradius, i.e., / m f p = r g = pc/{ZeB). 
Then 

Amfp = t/mfp Tg 5 (2) 

and the precursor diffusion length for species i is 

L l D oc r) mfp ji(A/Z)v? . (3) 

Here 7 , ( 17 ) is the Lorentz factor (velocity) for species i, 
c is the speed of light, and B is the background magnetic 
field in Gauss used to scale r g . If u, ~ c, then 7 oc p/A , 

L l D oc rj mlp (A/Z)(p/A) , (4) 

and the precursor diffusion length for different A/Z ions, 
at the same momentum per nucleon, scales as A/Z. 

Therefore, if ? 7 m f p is similar for all species, at the 
same p/A, high A/Z particles will diffuse further into 
the upstream region, ‘feel’ a larger effective compression 
ratio in the modified shock structure than low A/Z par¬ 
ticles, and gain a larger momentum boost in the next 
shock crossing. Since L l D increases with p/A , the modi¬ 
fied shock precu rsor produces a distincti ve concave spec¬ 
tral shape le.g.. lEllison fc Eichledfl~984fl with high A/Z 
particles having a harder spectrum at any p/A than low 
A/Z particles (see Fig. [13] below). It is important to 
note that equation J 2 | is presented only as a simple base¬ 
line scattering mode. Small-scale turbulence, as gener- 
ated by the Weibel ins t ability at relativ i stic shocks ( e.g., 
iNishikawa et al.1120061 : IPlotnikov et al.112011L 120131) , in 
general leads to a different A m f p oc p 2 dependence (e.g., 
iJokipiil fl971f) . In fact, the scattering process in rela¬ 
tivistic shocks is almost certain to be more complicat¬ 
ed than any simple power-law dependence for A m f p , and 
rj m fp may b e expected to have a momentum dependence 


as well (e.g.. 

Achterbers et al. 200l|Kirk & Revillc 2010fc 

Sironi et al.1 

20131: Lemoine et al.l 20141). Our current re- 


suits, with A m f p oc r g , display the essential physical ef¬ 
fects that come about from the development of an ex¬ 
tended shock precursor. As long as A m f p is an increasing 


function of p, and Fermi acceleration is efficient, the non¬ 
linear effects we describe should not depend qualitatively 
on the diffusion coefficient. Quantitatively, the momen¬ 
tum dependence of A m f p can strongly influence the max¬ 
imum CR energy a given shock can produce and this, 
in turn, will influence the photon production. We are 
currently generalizing our Monte Carlo simulation to al¬ 
low for a more complicated parameterization of A m f p over 
large energy and length scales. These futur e results w ill 
be co mpared against PIC simulations fe.g.. iSironi et all 
I2013H . The results presented here provide a benchmark 
for future comparisons. 

While the A/Z effect we describe can enhance the 
injection and acceleration of heavy ions compared to pro¬ 
tons, for electrons with A/Z ~ 5.45xl0 -4 it acts strongly 
in the opposite fashion. If only this kinematic effect is 
considered with equation electrons will be dramati¬ 
cally less efficiently injected and accelerated than protons 
in nonlinear Fermi acceleration. 

Electron injection was considered in a non- 
relativi stic Monte Carlo co de similar to the one we use 
here in iBaring et al.l (I1999T) . In that paper, in order to 
overcome the dramatic A/Z effect and allow electrons to 
be injected and accelerated with efficiencies large enough 
to be consistent with synchrotron and IC radiation ob¬ 
served in young supernova remnants (SNRs), the elec¬ 
tron A m f p was set equal to a constant below some momen¬ 
tum, i.e., changing rj m f p selectively for electrons. This 
modification gave low energy electrons a larger mean free 
path than equation ([2]) would produce and allowed them 
to diffuse far enough upstream to overcome the shock 
smoothing effects. It was argued in IBaring et al.1 (11999H 
that this simple modification was reasonably consistent 
with a n electron injection model developed by iLevinsonl 

PHD. 

Here, we adopt a different approach. We keep equa¬ 
tion ©, but transfer some fraction of the ram kinetic en¬ 
ergy from ions to electrons as the particles first cross the 
viscous subshock. We note that a “sharing” of energy be¬ 
tween ions and electrons is clearly seen in recent PI C rel a¬ 
tiv istic shock simulation s li.e.. ISironi h Spitkovskvll2011l; 
ISironi et al.l 12013 1, and IPlotnikov. Pelletier, k, Lemoinel 
(|2013fl give an analytical treatment of electromagnetic 
instabilities transferring energy to electrons in the pre¬ 
cursor of relativistic shocks. 

In the PIC simulations, electrons are heated in the 
precursor by interacting with turbulence generated main¬ 
ly by backstreaming protons and obtain near equiparti- 
tion with the protons before crossing the subshock. We 
mimic this effect by transferring a set fraction, /; on , of ion 
energy to electrons as particles first cross the subshock. 
While this simple energy transfer model is clearly an ap¬ 
proximation, we feel it affords a straightforward way of 
using plasma physics information obtained from compu¬ 
tationally intensive PIC simulations in a calculation that 
can model particle acceleration and photon emission con¬ 
sistent with the production of high-energy cosmic rays in 
relativistic shocksQ 

7 Other descriptions of electron heating in relativistic plasmas 















































Particle Acceleration in Relativistic shocks 


3 


We recognize, of course, that collisionless shock for¬ 
mation, and particle injection and acceleration, are de¬ 
termined by more than kinematics alone. The self gen- 
eratio n of magnetic turbulence is c ritica l to the process 
(e.g., iBednarz fe Ostrowskil [19981 iLcmoine fe Pelletier! 
[20Q3), particularly for relativistic shocks, and, as just 
mentioned, energy can be transferred between electrons 
and protons by wave-particle interactions that are es¬ 
sentially independent of the kinematics. Nevertheless, 
if Fermi acceleration is efficient, basic momentum and 
energy conservation demands that kinematics be taken 
into account and the shock precursor must be modified 
by the backpressure of accelerated particles. 

In contrast to non-relativistic shocks, where shock 
acceleration can be directly tested against spacecraft ob - 
servations fe.g.. lEllison et all 1 1 9901 [Baring et ~all 1 1 9971) . 
Fermi acceleration in relativistic shocks is far less certain. 
Relativistic shocks cannot be directly observed with 
spacecraft, as can non-relativistic, heliospheric shocks, 
and the highly anisotropic particle distributions intrinsic 
to relativistic shocks make the self-generation of mag¬ 
netic turbulence far more difficult to describe analytical¬ 
ly. Despite this difficulty, i ntensive work continues (e.g: ., 
iLemoine fc Pelletleil 1201(1120111 iPlotnikov et al.1 12013 1. 
Furthermore, and again in contrast to non-relativistic 
shocks, the predictions for particle spectra and, there¬ 
fore, photon signatures from relativistic shocks are high¬ 
ly uncertain. Sites harboring relativistic shocks, such as 
GRBs, can often be successfully modeled with alterna¬ 
tive acceleration mechani sms (e.g., magnetic reconnec¬ 
tion in the case o f GRBs; iMcKinnev fc IJzdenskvll2012t 
iSironi et al.ll2015h . weakening the link between Fermi ac¬ 
celeration theory and observation. 

Despite this uncertainty, there is compelling 
evidence, pr i mar ily from PIC sim ul ations (e.g., 
Hoshino et al.1 J1992 iKatol l2007t ISiro ni fc Sni tkovsk 


2009 : Kcshct et a l. 2009; Nishikawa et ah! 1201 


I 


Sironi fc Spitkovskvl l2011l ~ lStockcm et al.l l2012lj . that 


relativistic shocks do accelerate electrons and protons 
beyond the initial kinematic boost from a single shock 
crossing. These simulations also highlight the role of 
the magnetization parameter, cr B , in Fermi acceleration, 
where 


a B = 


B 2 0 


ATTnorripC 


2 ’ 


(5) 


Bq i s the upstream fie ld, and no is the number density 
fsee lBvkov fc Treumannl 120111 . for an alternative defini¬ 
tion_of_CT^} ; __^ _ __ 

ISironi et al.l (j2013l ) (see also fHaugbolld 120111 1 show 
that perpendicular electron-ion shocks with Lorentz fac¬ 
tors 70 < 150 inject and accelerate electrons and ions 
efficiently when cr B < 3x 10 -5 . In these weakly mag¬ 
netized plasmas, the self-generated turbulent field dom¬ 
inates the uniform B and the shock obliquity ceases to 
be important, an assumption often made in shock accel- 


includc IGcdalin et al . (2003), who invest igate the effects of a cross¬ 
shock potential, and I Kumar et al.l (120150 , who explore the electron 
heating behavior of Weibel-induced current filaments in a 2D PIC 
simulation. 


eration studies and one we make here with our plane- 
parallel shock assumption (c.f., Fig. E2). Quantifying this 
lack of dependence on the obliquity is particularly im¬ 
portant since ultra-relativistic shocks are essentially al¬ 
ways highly oblique—the R-field component along the 
shock face can be highly compressed—and strongly mag¬ 
netized oblique shocks are less able to inject and acceler¬ 
ate particles. We note that for typical interstellar medi¬ 
um (ISM) conditions, Bq ~ 3/tG, no ~ 0.03cm -3 , and 
so a B ~ 10 -9 . This value is orders of m a gnitu de be¬ 
low the threshold reported by ISironi et al.l (120131) , sug¬ 
gesting that ultra-relativistic shocks propagating in the 
normal ISM may be able to inject and accelerate elec¬ 
trons and ions far more efficiently than previously be¬ 
lieved, regardless of the shock obliquity. We note that for 
non-relativistic shocks the plasma /3 = nok B To/(B q/8tt) 
is the relevant parameter rather than equation ([5]) and 
oblique geomet ry may be imp ortan t for typical ISM pa¬ 
rameters fe.g., lOrlando et alJl201lD . Here fee is Boltz¬ 
mann’s constant and To is the ambient unshocked tem¬ 
perature. 

Given that weakly magnetized relativistic shocks 
can inject particles, the maximum energy these par¬ 
ticles obtain in a given shock remains uncertain , al¬ 
though arguments presented by ISironi et al.l (I2013D sug¬ 
gest that the acceleration reaches a maximum, T max , 

where T max /( 7 om p c 2 ) ~ ctg 1 ^ 4 . While this can be a sub¬ 
stantial energy (for a B = 3xl0~ 5 , T max ~ 0.2TeV), it 
is well below what is often assumed in suggesting that 
relativistic sho cks may produce ultra-h igh-energy cos¬ 
mic rays (e.g.. iKeshet fe Waxmanl 120051) . Here we sim¬ 
ply parameterize the maximum particle energy by set¬ 
ting a maximum shock size with a free escape boundary 
(FEB)0 

Next we descri be the gene r alizat ion of the Monte 
Carlo code used in lEllison et al.l (I2013D to include the in¬ 
jection and acceleration of electrons simultaneously with 
ions. The nonlinear shock structure is calculated includ¬ 
ing thermal leakage injection and the backreaction from 
all species. Full particle spectra are determined at var¬ 
ious positions relative to the subshock, along with the 
distributions of particles that escape at upstream and 
downstream FEBs. For given values of the ambient den¬ 
sity, magnetic field, and background photon field, we cal- 
culate the synchrotron emis sion using equation (6.7a) in 
iRvbicki fc Lightmanl (I1979D . the invers e Com pton (IC) 
emission using equation (9) in Uonesl (I1968D . and the 
pion-decay emission usin g pa rameterizations give n by 
iKamae et al.l (|2006l 120071) and lKelner et al.1 (120091) . For 
the pio n-decay from He 2+ , we use the scaling relation 
given in lBaring et al.l (I1999T) . Once the emission is deter¬ 
mined in the local frame it is transformed to the observer 
(i.e., ISM) frame. 

We note that radiation losses for electrons are only 

8 A FEB is a position beyond which particles are assumed to de¬ 
couple from the shock. Any actual shock will be finite in extent and 
at some point high-energy particles will obtain diffusion lengths 
comparable to the shock size and stream freely away. Our FEBs 
model a finite shock size within our steady-state, plane-parallel 
approximation. 




























































































4 


Warren, Ellison, Bykov, Lee 


considered during the acceleration process—once accel¬ 
erated, the electrons radiate without accounting for fur¬ 
ther losses. This so-called ‘thin target’ approximation, 
where the radiation length is assumed to be larger than 
the region between the upstream and downstream FEBs, 
is adequate for the steady-state examples given h ere and 
can b e relaxed in models of evolving GRBs (see IWarrenl 

USD- 

2. MODEL 

We assume the basic Monte Ca rlo scattering model, 
as described in lEllison et al.l (120131) for protons, applies 
equally to electrons and heavy ions. If so, all that is need¬ 
ed to describe the injection and acceleration of electrons 
is a parameter describing the transfer of energy from ions 
to electrons, mimicking t he effect seen in P IC si mulations 
of relativistic shocks (see lSironi fe Spitkovskvll201lD . 

Full details of the particle scattering model, thermal 
leakage injection, and the method for obtaining a self- 
consistent shock precursor structure when Fermi accel- 
er ation and par ti cle es cape occur are given in S ection 2 
of lEllison et all (I2013T) . In lEllison et al.1 (120131) we al¬ 
so fully explain the caveats needed when applying the 
Monte Carlo model to relativistic shocks and show how 
the approximations required in the Monte Carlo model 
compare with previous Monte Carlo work and with more 
fundamental PIC simulations. 

We describe particle transport by assuming the mean 
free path, A m f P , is given by equation (J2]) with ? 7 m f p = 1, 
i.e., Bohm diffusion^ As mentioned above, while there is 
some theoretical and observational sup port for Boh m dif¬ 
fusion from X-ray afterglows (e.g.. iSagi fc Nakarll2012 :). 
most theoretical work suggests a stronger momentum de - 
pendence, i.e., Amfp oc p 2 (e.g.. iPlotnikov et al.l 1201 J) . 
We use equation © for convenience and as a baseline for 
future work which will assume more complicated (and, 
hopefully, more realistic) forms for A m f p . The particle is 
moved for a time St Ct c , where t c = A m f p /vi is the “col¬ 
lision time,” i.e., the average time measured in the local 
frame needed for the particle to accumulate deflections 
on the order of 90°. The second scattering parameter, 
N g , determines the “fineness” of scattering through an 
equation for the maximum deflection, S6 ma . x , a particle 
can experience in a pitch-angle interaction event after 
each St, i.e., 

56 

max — \J Vlit/(riroSpNg) . ( 6 ) 

Here, N g = T g /St is the number of gyro-time segments 
St dividing a gyro-period r g = 2irr g /vi, and we note 
that equation (O applies even if particles move rapidly 
between inertial frames, as is normally the case for rel¬ 
ativistic shocks. In each scattering event, the scattering 

9 We note that since all lengths in the steady-state code are 
scaled with r g o = Cm f P m P uo c/(e Bq) , where «o is the shock speed, 
our results are independent of r) m f p except for the absolute normal¬ 
ization of the particle spectra, i.e., the number of particles within 
a physical region scales as Cmfp- In ou r plane-parallel approxima¬ 
tion, the field throughout the shock retains the far upstream value 

B 0 . 


is assumed to be isotropic and elastic in the local plasma 
frame. 

Large values of N g imply fine scattering while small 
values imply that the particle momentum will suffer rel¬ 
atively large deviations i n direct ion i n each pitc h-angle 
scattering event (e.g., iSummerlin fc Baring! 120128) . In 
all examples here, N g is large enough to saturate and 
produce fine-scattering results that do not change sub¬ 
stantially as N g is increased further^ It is important 
to note that while we can simply parameterize particle 
transport with the assumption given by equation 0 and 
the two parameters A m f p and N g , this assumes that mag¬ 
netic fluctuations with correlation lengths on the order 
of L c = 2Trri m i p r g /Ng exist with sufficient power to pro¬ 
duce this scattering throughout the shock. We make 
no claim of self-consistently determining the magnetic 
turbulence needed to produce the diffusion implied by 
equation 00 

All of the above is applied equally to electrons and 
ions and this parameterization yields injection rates and 
acceleration efficiencies that depend on the AjZ ratio. 
We note that particle injection from shock-heated ther¬ 
mal particles (our thermal leakage injection model) oc¬ 
curs directly from the above assumptions. Regardless of 
the shock speed, virtually all cold, thermal, unshocked 
particles obtain a velocity greater than the downstream 
bulk flow speed when they first scatter in the downstream 
region and the flow becomes subsonic. Depending on 
their angular distribution, these particles have a finite 
probability of scattering back upstream where they will 
be further accelerated. In the Monte Carlo code, the 
injection efficiency is determined stochastically as some 
particles manage to scatter back upstream with an an¬ 
gular distribution determined solely from the scattering 
model with no additional assumptions or parameters oth¬ 
er than that the subshock is assumed to be transparent 
for back-scattering particles. 

To this we add one additional parameter, /i on , i.e., 
the fraction of far upstream ion ram kinetic energy 
transferred to electrons. The ion ram kinetic energy 
is defined as (70 — ljm^c 2 , where rrii is the ion mass, 
70 = [1 — (ito/c ) 2 ] -1 / 2 is the far upstream shock Lorentz 
factor, and uq is the far upstream (i.e., unmodified) 
shock speed. When an ion crosses the subshock from 
upstream to downstream for the first time, energy equal 
to /ion (70 — 1 )rriiC 2 is removed from it. The total ion 
energy transferred is /i on )Ci( 7 o — l)A/mfC 2 , where the 
sum is over the ion species and N t is the number of i- 

10 We note that all of our fully relativistic, unmodified shock 
examples yield the well known power law dN /dp oc p~ 2 23 (e.g., 
Bednarz fc Ostrowskil I i 99 Si : iKirk et al.l 12000 :: Kcshct fe Waxmanl 
200oTi (see Fig. [3] below). Harder spectra can be obtained 
with large-angle scat t ering, i.e., small values of N g (e.g., 
lEllison fe Double! |2004|; [S ummerlin 7 B aring) 1201211 , or other ad 
hoc assumptions (e.g.. ISchlickeiserll2015ll but care must be taken 
to account for NL effects when the power law index is less than 2 
(see figure 13 in lEllison et al.1120131) . 

11 We note that a nonlinear, non-relativistic Monte Carlo code 
which does include th e self-generation of mag n etic t u rbulence has 
been developed (e.g., Vladimirov et al. 2008. 2 0091 : IBvkov et al.l 
120141) . 
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Fig. 1.— Schematic, not-to-scale representation of a relativistic 
shock embedded in a conical jet, propagating with a Lorentz factor 
70 into material at rest. The region of interest is delimited by LypS 
and L Dw s; outside of this region, we assume the particles decouple 
from the shock and we ignore any emission resulting from them. 
The jet opening half-angle is djj , and the shock has propagated a 
distance i? s hk e f in the explosion (or ISM) frame. We only consider 
emission directed along the jet to an observer at Earth. 



Fig. 2.— This detail focuses on the volume between I/UpS and 
LdwS- The upper panel shows the shock-frame plasma velocity 
profile. The lower panel shows the ISM-frame velocity at select¬ 
ed points in the shock structure, varying from 0 for x < LupS to 
u o just upstream from the subshock at x = 0. Note that we as¬ 
sume the shock is locally plane and that Lims and TdwS are small 
compared to iJ s hk,ef as illustrated in Fig. [l] 

ions injected far upstream. This ion energy is divided 
equally among electrons and added to their energy as 
they cross the subshock into the downstream region for 
the first time. The parameter /j on can increase electron 
injection substantially as we show below. 

3. RESULTS 

We approximate the geometry of a relativistic after¬ 
glow shock moving in a jet as shown in Fig. [Q A detail 
of this situation is shown in Fig. [2j where we show the 
additional approximations for this preliminary work that 
the shock is locally plane, the distances to the upstream 
or downstream FEBs (TupS and Tdws) measured from 
the subshock are small enough so the jet cone is approx¬ 
imately a cylinder in the region surrounding the shock, 
and the diameter of the “cylinder” is large enough so 
particle escape out the sides of the cone is negligible. 
The cone material outside the upstream FEB is assumed 
to be stationary, i.e., it is in the local ISM or explosion 
frame. Outside of the region between TupS and TdwS 
we assume all accelerated particles have decoupled from 
the plasma and we ignore any emission they might pro¬ 
duce. A more realistic GRB afterglo w model follow ing 
the evolution of a jet shock is given in IWarreiil (120151) . 

With these approximations we calculate the shock 



L°g 10 P [m p c] 


Fig. 3.— Protons (black curves) and electrons (red curves) from 
UM shocks with different f[ on as indicated. These spectra, multi¬ 
plied by p 2-23 , are calculated downstream from the shock, in the 
shock frame, and have arbitrary overall normalization although 
the relative normalization between electrons and protons is abso¬ 
lute. An upstream FEB of Lu p s = —10 4 r^o was used with no 
downstream FEB, i.e., a probaility-of-return calculation was used 
to simulate an infinite downstream region. In Table fTI the /i on = 0 
case is Model A and the f[ on = 0.15 case is Model B. 



L °gio p [%c] 

Fig. 4.— Proton and electron spectra (as labeled and multiplied 
by p 223 ) with /i on = 0.15. These spectra are measured in the local 
plasma frame and are normalized to the total number of particles 
in a given region, as indicated in Fig. [2] The shock acceleration is 
limited by an upstream FEB at LupS = — 1 000 r g o : and a down¬ 
stream FEB at LdwS — +1000 r g o- 


structure and particle spectra for a given set of param¬ 
eters, as listed in Table [l] First we consider unmodi¬ 
fied (UM) shocks where the backreaction of the accel¬ 
erated particles on the shock structure is ignored. We 
note an essential difference between unmodified shocks 
and test-particle (TP) ones. A TP shock is one where 
the injection and acceleration efficiencies are low enough 
so the backpressure from accelerated particles can be ig- 
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Fig. 5.— The top panel shows the shock structure for the un¬ 
modified case (dashed red curve, Model D) and the nonlinear 
case (solid black and dot-dashed blue curves, Model E). The solid 
(black) and dashed (red) curves in the top panel are the flow speed 
u{x)/uq while the dot-dashed (blue) curve is 7 ( x ) 3 (.x) / ( 70 3 0 ) , 
where x(x)8(x) scales as 1/density for the nonlinear shock. The 
middle and bottom panels show the momentum and energy fluxes, 
respectively, normalized to far upstream values. Note that three 
particle species, protons, He 2 +, and electrons, are included in de¬ 
termining the self-consistent shock structure. 


nored. In the TP case, to the limit of total energy placed 
in accelerated particles, momentum and energy can be 
conserved without modifying the shock structure. Our 
thermal leakage injection model is efficient enough so TP 
shocks are never produced for the parameters we use 
here. In our UM examples, particles are injected and 
accelerated efficiently but the effect of shock accelerated 
particles on the shock structure is ignored. This provides 
a direct comparison to nonlinear (NL) shocks, where mo¬ 
mentum and energy are conserved. Our main point in 
this paper is that, if the acceleration is efficient, the shock 
must be modified by the accelerated particles. 

While we do not consider TP shocks explicitly, the 
superthermal particle fluxes and photon emission from 
our UM examples can be simply re-scaled to a TP re¬ 
sult. For example, if a TP result is defined as one where 
superthermal particles contain < 1 % of the total energy 
flux, so energy flux will be conserved to within 1 %, our 
UM superthermal particle and photon emission needs to 
be reduced by a factor /Vtp > 100 x [F en (x > 0 )/F e ° n ], 
where F e ° n is the far upstream energy flux and F ea (x > 0) 
is the UM downstream energy flux. For the UM 70 = 10 
examples we show below, F en (x > 0)/F e ° n ~ 2. 

3.1. Unmodifed (UM) Examples 


In Fig. [3]we show number spectra, dN/dp, calculated 
just downstream from an unmodified shock with 70 = 10 . 
The solid and dotted black curves are protons and the 
red curves are electrons. The dotted curves were calcu¬ 
lated for /i on = 0 , while the solid curves were calculated 
with /i on = 0.15, i.e., 15% of the proton ram kinetic ener¬ 
gy is transferred to electrons as particles first cross x = 0 
headed downstream J3 Since all particles injected far up¬ 
stream are non-relativistic, the transformation from up¬ 
stream to downstream frames strongly favors more mas¬ 
sive particles in the first shock crossing. This results in 
the strong depression of electrons relative to H + . Once 
all particles become relativistic they are treated equal¬ 
ly and obtain similar power laws, i.e., the UM shocks 
in Fig. [3] show the canonical dN/dp oc p~ 223 power law 
above the “thermal” peak and below the high momen¬ 
tum cutoff. For the protons, the cutoff at ~ 10 7 m p c is 
produced by an upstream FEB at LupS = -lx 10 4 r s o, 
where r g 0 = r] m { p m p uoc/(eBo). The electrons cut off 
at a lower momentum (~ 2 x 10 6 m p c ) due to radiation 
losses. Without radiation losses, the electrons would ob¬ 
tain the same p max as protons since p max scales as Z. 
The normalization of the electron spectra shows the dra¬ 
matic effect of fi on : the e/p ratio is increased by nearly 
three orders of magnitude with /i on = 0.15. The proton 
normalization is only slightly influenced by /; on . 

To save computation time, the models in Fig. [3]used 
a probability-of-return calculation instead of a down¬ 
stream FEB. This mimics an infinite downstream re- 
gion and a ll ows ra pid acceleration to high energies (see 
lEllison ct a 1.1 1 199(1 for a discussion of the probability of 
return calculation). In Fig. 2] we show spectra from an 
UM shock measured in the local plasma frame with both 
upstream and downstream FEBs. While both FEBs are 
present, typically one (the shorter, measured in diffu¬ 
sion lengths) will determine the maximum momentum 
Pmax- For Fig. IH both L UpS and L D wS are 1 x 10 3 r g0 
from the subshock, but LdwS determines p m ax since it 
is much more difficult for particles to stream away from 
the shock in the upstream region with 70 = 10 . 

The FEBs also determine the total number of parti¬ 
cles accelerated (as indicated in Fig. m and the dN/dp 
spectra in Fig.[4]are normalized to the total particle num¬ 
ber in a region surrounding an observation position x 
(the position is indicated in Fig. 0 • While the particular 
number of accelerated particles in these examples is arbi- 
trarvFl it is important to note that the Monte Carlo code 
determines the absolute number of accelerated particles, 
and subsequent radiation, for any given set of environ¬ 
mental and shock parameters, including /i on - As long as 
70 3> 1 , the large majority of accelerated particles will be 
in the downstream region since it is difficult for particles 
to stream upstream into the shock precursor. This is re- 

12 We refer to the sharp d rop in u{x) that occurs at; x — 0 as 
the subshock (see Figs.[5land ll0ll . For an unmodified shock, there 
is no distinction between the shock and subshock. 

13 The normalization of the spectra in Fig. [4] (and Figs. [6] and 
1111 below) depends on the opening angle and length of the jet, 
t he ambient number density, and the positions of the FEBs (see 
lWarrcnll2015l . for a full discussion). 
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Fig. 6.— Downstream, local plasma frame (LPF) spectra for the 
unmodified shock shown in Fig. [5] (top panel, Model D) and the 
nonlinear shock shown in Fig. [5] (bottom panel, Model E). Note 
the pronounced “superthermal” tail on the electron distribution. 


fleeted in the higher normalizations of the black curves, 
measured at x = 0 + (i.e., downstream), compared to the 
red or blue curves, measured in the shock precursor at 
x = — 1 7'go and x = —120 r g o, respectively. 

Besides the downstream FEB and normalization, the 
spectra in Fig. 0] differ from those in Fig. [3] in that they 
are shown in the local plasma frame and they include 
spectra calculated upstream from the subshock (red and 
blue curves), as well as downstream (black curves). The 
effect of the Lorentz transformation from the shock 
to the local plasma frame is clearly indicated by the 
high-momentum upstream spectra (red and blue curves) 
which extend to higher momentum than the downstream 
spectra. Relative to the subshock, the upstream plasma 
frame moves with 70 = 10 , while the downstream plasma 
frame moves with 72 = [1 — (it 2 /c ) 2 ] _1 / 2 ~ 1.06. Here 112 
is the downstream bulk plasma speed measured in the 
shock frame. 


3.2. Nonlinear Examples , 70 = 10 

In Fig. Owe show the structure of a shock where the 
backpressure from accelerated protons, He 2+ , and elec¬ 
trons is taken into account. For this example, f lon = 0.1, 
and the electron pressure contributes to the determina¬ 
tion of the self-consistent shock structure. Here AupS = 
—300 r g o and LdwS = +1000 r g o- In the top panels, the 
solid (black) curve is the bulk flow speed, u(x)/uo, and 
the dot-dashed (blue) curve is r y(x)/3(x)/('yofio), where 
P(x) = u(x)/c , 7 (x) = [1 - /3{x) 2 )~ 1/2 , and /3 0 = u 0 /c. 



Log 10 MeV 

Fig. 7. — Photon emission for the UM (Model D) and NL (Model 
E) shocks shown in Figs. f5l and l6l The dashed (black) curves are 
the total emission transformed to the ISM frame for an observer 
at -Dobs = -c ^Mpc Mpc, at an angle within I /70 from the shock 
normal. The lower synchrotron and pion-decay curves in each panel 
show emission calculated in the local plasma frame and summed 
over the shock from Lu p s = —300r^o to I/DwS — +1000 r g Q. The 
dashed (red) pion-decay curves are the LPF emission from He 2-1- . 

The dashed (red) curve is u{x)/u o for the UM case. The 
lower panels show the momentum and energy fluxes for 
the NL case (solid black curves), as well as for the UM 
case (dashed red curves). All curves are normalized to far 
upstream values. Without shock smoothing, the down¬ 
stream momentum and energy fluxes are nearly a factor 
of two out of conservation. 

We note that in non-relativistic and trans-rela- 
tivistic shocks, the smoothing required for momen¬ 
tum and energy conservation is accompanied by an 
increase in the overall shock compression ratio above 
the Rankine- Hu goniot v alue, i.e., i?tot > -Rrh (e.g., 
[Berezhko fo Ellisonll19991 ). The overall compression ratio 
is defined as I?tot = uq/u 2 and the subshock compres¬ 
sion ratio is R s u b = u\/u 2 - Here, 112 is the bulk plasma 
speed downstream from the subshock and u\ is the plas¬ 
ma speed just upstream of the viscous subshock, both 
measured in the shock rest frame F^l 

14 As seen in Fig. [5] the definition of ui is imprecise because the 
Monte Carlo solution allows for a smooth decrease in the precur¬ 
sor speed into the downstream region. Part of the increase in ittot 
is due to escaping particles acting to soften the equation of state 
for the plasma that remains coupled to the shock system. The re- 
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The top panel in Fig. 0 shows downstream proton, 
He 2+ , and electron spectra for the UM shock (dashed 
curves in Fig. 0. The bottom panel shows these spectra 
from the nonlinear shock. All parameters are the same 
for these two cases—the only difference is that momen¬ 
tum and energy are conserved in the NL case. 

The effects of the smooth shock structure are clear¬ 
ly evident. In the NL case, the downstream spectra are 
noticeably curved and less intense, as is necessary to con¬ 
serve energy and momentum. Significantly, the electrons 
are more modified than the protons or He 2+ and the e/p 
ratio, in the quasi-power law portion of the spectra, drops 
by more than an order of magnitude. With /; on = 0, e/p 
would have dropped by several more orders. This dif¬ 
ference is a direct result of our scattering assumption, 
i.e., equation ©. The electrons, with their small A/Z , 
feel the effects of the smooth shock more acutely than 
the heavier ions and are less efficiently injected and ac¬ 
celerated until they reach p > 10 m p c. In sources where 
electrons are assumed to have a large fraction of the ener¬ 
gy budget, the efficiency for accelerating electrons cannot 
be orders of magnitude less than it is for protons. There¬ 
fore energy must be transferred from heavy particles to 
electrons with a reasonable e fficiency if Fer mi accelera¬ 
tion is to be impo rtant (e.s.. ISironi fe Spitkovskvll20lH 
iSironi et al.ll2013l) . 

The heavier He 2+ , with A/Z = 2, is accelerated 
more efficiently than protons; the He +2 /p ratio above 
~ 100 TOpC goes from He +2 /p < 1 in the UM shock to 
He +2 /p > 1 in the NL shock (see Fig. [T51) . This is par¬ 
ticularly significant at the high momentum cutoff. The 
A/Z effect is discussed in more detail in Section I5T)1 For 
Model E in Fig. 0 the fraction of total ram kinetic ener¬ 
gy placed in particles of 100 MeV or greater, is £h =0.60, 
£He = 0.30, and e e i = 0.10, for protons, He 2+ , and elec¬ 
trons, respectively (see Table 0 where it is noted that 
these fractions are measured in the shock frame.) 

The effects of shock smoothing also show up in the 
broadening and shift to lower momentum of the “ther¬ 
mal” peaks, as seen in the bottom panel of Fig. 0 Of 
particular interest is the pronounced “superthermal” tail 
the electrons obtain in the NL case. If synchrotron-self- 
absorption (SSA) is unimportant, this can produce a no¬ 
table effect in the synchrotron emission, as we discuss 
next. 

3.3. Photon Emission, 70 = 10 

In Fig. 0we show the photon emission produced by 
the shocks described in Figs. 0 and 0 As in Fig. 0 
the top panel is for the unmodified shock and the bot¬ 
tom panel is for the nonlinear shock- all other parame- 

mainder of the increase is caused by the ratio of specific heats for 
particles downstream of the shock: as the shock speed increases, 
the particles receive more energy from the first shock crossing, the 
average particle approaches relativistic energies, and the ratio of 
specific heats drops fro m 5/3 —> 4/3. For fully relativistic shocks, 
however, Rtot — Rrh JDouble et al.|[2004J) : there is minimal par¬ 
ticle escape, and the downstream ratio of specific heats is already 
~ 4/3 and so cannot decrease further. The compression ratio in 
Fig. [5] reflects this invariance, as it is the same in the unmodified 
and nonlinear cases. 


ters are the same. The curves labeled synchrotron and 
pion-decay show isotropic emission calculated in the lo¬ 
cal plasma frame (LPF) summed over the regions be¬ 
tween the upstream and downstream FEBs (i.e., between 
x = —300 r g o and x = +1000 r g o) as indicated in Fig. 0 
These are the fluxes that would be observed at a distance 
dyip C Mpc if no Lorentz transformations were required. 
Here, dM pc is the distance in Mpc. Of course the parti¬ 
cle distributions will not be isotropic in the ISM frame, 
and Lorentz transformations are required. The dashed 
(black) curves show the total emission from synchrotron, 
IC, and pion-decay transformed to the ISM frame and 
seen by an observer at D a b s = — drupe Mpc within an 
angle 1 /70 from the shock normal. 

The shock simulation is done in the shock rest frame. 
To obtain the local plasma frame synchrotron and pion- 
decay spectra we first transform the particle spectra to 
the LPF taking into account the anisotropies introduced 
by the relativistic flow. We then calculate the photon 
emission in the LPF assuming it is produced isotropical¬ 
ly. This is reasonable for the synchrotron emission since 
we implicitly assume the background magnetic field is 
highly turbulent, and the synchrotron photons should 
be produced isotropically as the electrons spiral in the 
turbulent field. It is also a good approximation for the 
pion-decay emission since the protons and He 2+ ions in¬ 
teract with the local plasma and produce pions that can 
become isotropic in the LPF before emitting a 7 -ray. 

From this isotropic emission, we obtain the flux in 
the ISM frame ahead of the shock by employing the stan- 
dard Doppler shift an d Lorentz transformations (e.g., 
iKumar fc Zhand I2014T) . For an observer at D a b s = 
—dMpcMpc, within an angle 0 o bs = I /70 from the jet 
direction, these give a boost to the energy flux oc £> 4 , 
where V = j(x) [1 + /3(x) cos0 o bs] is the Doppler factor 
at position x relative to the subshock. One factor of V 
comes from the Doppler shift, one from time dilation, 
and two from the relativistic beaming. Here y(x) and 
/3(x) are measured relative to the ISM frame and will be 
different for each region between the upstream and down¬ 
stream FEBs as shown in Fig. 0 Far upstream T> = 1, 
while downstream for 70 = 10 and Rtot — 3.02, V ~ 14. 

Since the cosmic microwave background (CMB) pho¬ 
tons are nearly isotropic in the ISM frame, we calculate 
the IC emission in a different fashion{0 We first trans¬ 
form the electron distribution into the ISM frame keeping 
the two-dimensional anisotropy inherent in our plane- 
parallel shock simulation. We then calculate the emis¬ 
sion directed toward the observer at D a b s = —^MpcMpc 
using equation (9) in lJonesI (11968H . The assumption here 
is that the relativistic electrons produce strongly beamed 
emission, so only those electrons directed toward the ob¬ 
server contribute to the observed flux. In Fig. 0 the 
dashed (black) curves contain the full observed flux from 
emission produced over the modified shock structure and 
transformed to the ISM frame for an observer within 

15 We only consider CMB photons here for simplicity—the tech¬ 
niques we present can be generalized to include IC emission from 
other photon fields including synchrotron self-Compton photons if 
the jet parameters warrant it. 
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Fig. 8 .— Nonlinear downstream LPF spectra for Model E (/i on — 
0.1) and Model F (/i on = 0.4). In both panels the solid (black) 
curves are protons, the dashed (red) curves are He 2- * - , and the 
dotted (blue) curves are electrons. 


I /70 from the jet axis (i.e., toward -\-x in Fig. [2j). The 
IC emission is part of this sum. 

The effects from shock smoothing on the particle dis¬ 
tributions (Fig. [ 6 ]) produce corresponding changes in the 
photon emission. Since electrons are suppressed more 
than ions in the NL shock, the synchrotron and IC emis¬ 
sion drops more than pion-decay between the UM and 
NL cases. In contrast, the pion-decay from He 2+ (dashed 
red pion-decay curve) is increased relative to that from 
protons by the A/Z effect. The “thermal” peaks near 
1 rripC for electrons, and near 10 m p c for protons and 
He 2+ , show up as clear peaks in the synchrotron and 
pion-decay emission at E ~ 10 -11 and 10 3 MeV, respec¬ 
tively. This is particularly significant for the synchrotron 
emission near 10 ~ 9 MeV. The NL curvature in the par¬ 
ticle spectra shows clearly in the individual components 
and remains strongly evident in the summed flux (bot¬ 
tom panel of Fig. 0. 


3.4. Nonlinear, 70 = 10, Variation of fi on 

While we have us ed f; nn = 0 .1 in ou r NL Model E, 
the PIC simulations of lSironi et al.l (120131 ) show examples 
where ~ 40% of the energy in ac celerated particles ends 
up in electrons (see figure 11 in iSironi et al.ll20l3l) . In 
Fig. [ 8 ] we compare particle spectra, and in Fig. 0 the 
total observed energy flux, for /j on = 0.1 (Model E) and 
/ion = 0.4 (Model F). The flux between ~ 10 -9 MeV and 
~ 1 GeV is ~ 100 times greater for the /i on = 0.4 case, 



Log 10 MeV 


Fig. 9.— Total observed energy flux for NL Models E (/i on = 0.1) 
and F (/ ion = 0.4). 
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Fig. 10.— All curves are as in Fig. [5] for 70 = 1.5. The 
dashed (red) curves are the UM (Model G) case while the solid 
(black) and dot-dashed (blue) curves are the NL (Model H) pro¬ 
files. The upstream and downstream FEBs are LupS = — 300 r g o 
and LdwS = +1000 r g o, as in Fig. [5] but these distances differ in 
absolute units since uo varies between 70 = 10 and 70 = 1.5. 


with a much smaller decrease in the GeV-TeV emission. 
The /i on = 0.4 example shows the curved spectral shape 
that results from the NL shock smoothing but it is less 
pronounced between 10 -9 and 1 MeV than for /; on = 0.1. 
In the GeV-TeV range, the curvature is slightly greater 
than for / ion = 0.4. 

With /ion = 0.4, the energy distribution above 
100MeV in protons, He 2+ , and electrons is en = 0.36, 
eHe = 0.20, and e e i = 0.44, respectively; nearly 50% 
of the ram kinetic energy goes into 100 MeV or greater 
electrons, as measured in the shock frame. The large dif- 
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Fig. 11.— Except for the dotted (black) curves labeled ‘70 — 10 
p’s’ all spectra are downstream, local plasma frame (LPF) spec¬ 
tra for the unmodified 70 = 1.5 shock (top panel, Model G) and 
the non linear 70 = 1.5 shock (bottom panel, Model H) shown in 
Fig. 1101 The curves labeled ‘70 = 10 p’s’ are identical to those in 

Fig- El 


ference in electron normalization between the /i on = 0.1 
and 0.4 cases comes about mainly from the NL shock 
smoothing effects the electrons with /; on = 0.1 receive. 
When / ion = 0.4, the lowest energy downstream elec¬ 
trons that cross back upstream have a long enough dif¬ 
fusion length so they feel a large effective compression 
ratio. Nonlinear effects reduce the electron acceleration 
more for f- lon = 0.1 than for /i on = 0.4. 


3.5. Trans-relativistic, 70 = 1.5 

As described in lEllison et all (120131) . the Monte Car¬ 
lo simulation smoothly treats non-relativistic to ultra- 
relativistic shocks. In Fig. [TO] we show the profile of 
a trans-relativistic 70 = 1.5 shock for comparison with 
Fig. [5] Apart from 70 , all input parameters are the same 
for the 70 = 1.5 and 70 = 10 cases. Since it is easier 
for accelerated particles to diffuse upstream against the 
inflowing plasma with 70 = 1.5 than against 70 = 10, 
the NL shock precursor is much more extended than it 
is for 70 = 10. The bulk flow speed, u{x ), is noticeably 
modified out to x = Tups = —300 r s o (note the split 
log-linear x-axis in Fig. fTtil) . 

As in the 70 = 10 case, the momentum and ener¬ 
gy fluxes are not conserved in the UM shock but are 
within a few percent of the far upstream values once 
the shock structure is modified by the CR backpressure. 
For 70 = 1.5, the overall compression ratio must al¬ 
so be increased to conserve momentum and energy and 



Log 10 MeV 

Fig. 12.— Photon emission for the U M ( Model G) and NL (Mod¬ 
el H) shocks shown in Figs. 1101 and 1 11 1 Except for the dotted 
(black) curves labeled ‘Total, ISM 70 = 10’, all curves are for the 
70 = 1.5 shock and are in the same format as in Fig. [7] The dashed 
(black) curves are the total emission transformed to the ISM frame 
for an observer at x = —c^Mpc Mpc, at an angle within I /70 from 
the shock normal. The dot-dashed (blue) curve in the lower panel 
is the total ISM frame 70 = 1.5 emission multiplied by 10 5 . 

we find Rtot = 3.9 ± 0.4, where the uncertainty comes 
from statistics and errors inherent in the Monte Car¬ 
lo smoo thing algorith m. This result is similar to that 
given in lEllison et all (I2013D except here we have includ¬ 
ed He 2+ and electrons in determining the self-consistent 
shock structure, and we use a downstream FEB as well 
as an upstream one. 

In Fig.[TT]we show the particle spectra for 70 = 1.5 in 
the same format as Fig. [HI except we have added the 70 = 
10 proton spectra for comparison. The 70 = 1.5 spectra 
are harder than those for 70 = 10 mainly because Rtot 
is larger. A comparison of the e/p ratio at ~ 300 m v c 
for the NL models in Fig. [ 6 ] ( e/p ~ 1/300) and Fig. [fl] 
{e/p ~ 1/6), shows that the 70 = 1.5 shock is much more 
effective in injecting and accelerating electrons than the 
70 = 10 shock. However, the maximum momentum is 
noticeably lower for 70 = 1.5 and the “thermal” peak 
is also lower since, for downstream spectra, it occurs at 
~ 70 m p c for protons in the UM shock. The “thermal” 
peak is at a noticeably lower momentum in the NL case. 

Another important difference is the normalization. 
The 70 = 10 proton distributions are about a factor of 
100 above the 70 = 1.5 spectra in the quasi-power law 
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Fig. 13.— Models I, D, and G show proton and He 2 "*" spectra 
for UM shocks, while models J, E, and H are the corresponding 
spectra for the NL shocks. In all cases, the He 2+ spectra are mul¬ 
tiplied by 10 to adjust for the ambient number density. When 
plotted in p/A units the UM spectra are identical except for statis¬ 
tics and the maximum momentum cutoff. The NL shocks show a 
clear A/Z enhancement in the He 2 +/p ratio. 

region. This comes about for two primary reasons. The 
7 o = 10 shock has a considerably higher downstream 
density since ni = 7o/3o?io/( 72 /^ 2 ) and, with r g 0 be¬ 
ing larger, the number of particles in the distribution 
is larger. Another cause is that the spectra are mul¬ 
tiplied by p 223 rather than by p. The area under a 
dN/d\og{p) curve would be the total number of parti¬ 
cles. But dN/dlog(p) = p ■ dN/dp , so the extra 1.23 
powers of p merely enhance the perceived importance 
of particles at higher momentum. Since the 70 = 10 
shock contains particles at higher momenta, it also gets 
plotted higher on the p 2 23 dN/dp plot. The efficiencies 
for producing 100 MeV or greater energy particles are 
£h = 0.52, €He = 0.27, and e e i = 0.08 for Model H; i.e., 
the 70 = 1.5 shock puts ~ 8 % of the shock energy into 
energetic electrons. 


In Fig. 0 we show the photon emission for 70 = 1.5 
in the same format as Fig. 0 with the addition of the 
total ISM frame emission for the 70 = 10 cases (dotted, 
black curves) added to the 70 = 1.5 cases (dashed, black 
curves) for comparison. The total emission for an observ¬ 
er at Hobs = — (iMpc Mpc is dramatically different for the 
two Lorentz factors. To highlight this in the lower pan¬ 
el, we have plotted the total ISM emission for 70 = 1.5 
multiplied by 10° (dot-dashed, purple curve). The broad¬ 
band spectral shapes are very different between 70 = 10 
and 1.5, and at ~ IkeV the normalization differs by 
~ 10 5 . 

Note that the NL 70 = 1.5 spectrum is harder than 
70 = 10 in the GeV-TeV range but cuts off at a lower en¬ 
ergy. The lower cutoff energy shows up dramatically in 
the synchrotron emission. For 70 = 1.5 the synchrotron 
peak is around 1 keV, typical of SNR observations, while 
for 70 = 10, the peak is around 1 MeV. At radio emitting 
energies, the synchrotron spectra are very different be¬ 
cause of the emission produced by the downstream ther¬ 
mal electrons. Of course we have not considered SSA 
here and this process will produce a low-energy cutoff in 
the synchrotron emission which may mask the emission 
from the thermal electrons. Since the frequencies of the 
thermal peak and SSA cutoff depend on B as vt\i oc B 
and v a oc I? 1 / 5 , respectively, a stronger magnetic field 
makes it more likely that the thermal peak will be re¬ 
solved. We have also not considered 7 -ray absorption 
between the GRB and Earth. 

3.6. A/Z Enhancement of Heavy Ions 

The Monte Carlo code assumes that all scatterings 
are elastic in the local plasma frame. This implies that 
an insignificant fraction of the particle energy is trans¬ 
ferred to magnetic turbulence in the wave generation 
process. With elastic scattering, 7,17 oc Pi/A remains 
constant in a scattering event, where 7 , and 17 are local 
frame values. The energy gain particles receive on cross¬ 
ing the shock, which is determined by a Lorentz trans¬ 
formation between the two frames, also scales as pt/A. 
In our plane-parallel approximation,the probability that 
particles make a set number of shock crossings also de¬ 
pends only on 7 , 17 . Thus an UM shock, with the Monte 
Carlo assumptions, will treat all particles identically in 
momentum per nucleon, including the thermal leakage 
injection. 

An exception to this occurs if the acceleration is lim¬ 
ited by a boundary at a fixed distance, as we assume 
here, or by a maximum acceleration time. Since we as¬ 
sume equation 0 , diffusion length and acceleration time 
both scale as ( [A/Z)(pi/A ) for Vi ~ c. Thus, apart from 
the normalization set by input parameters and the maxi¬ 
mum momentum cutoff, all species should have identical 
spectra when plotted against p/A. With a fixed FEB, 
high A/Z particles will turn over at a lower p/A than 
low A/Z particles. 

For models I, D, and G in Fig. [131 we show down¬ 
stream, shock frame, proton and He 2+ spectra for UM 
shocks plotted in p/A units. Except for statistical varia¬ 
tions, a factor of 10 normalization since n a = 0 . 1 n p , and 
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the high momentum cutoff, the proton and He 2+ spec¬ 
tra are identical. Electrons are not plotted but would 
show the same effect. In the corresponding NL models 
J, E, and H, a clear enhancement of He 2+ , produced 
solely because the shock structure in smoothed by the 
backpressure of accelerated particles, is seen. For the 
7 o = 30 Model J, cr = 0.66 and £R e = 0.34. 

The factor of two enhancement in the He 2+ /p ratio 
seen in the NL models in Fig. [13] should be large enough 
to see clearly in PIC simulations. Since this enhance¬ 
ment is a prediction that stems directly from important 
assumptions of efficient Fermi acceleration and thermal 
leakage injection, adding helium to PIC simulations can 
test these assumptions. If the acceleration is efficient, 
and the A/Z effect is not seen, it implies that one or 
more of the following may be happening, (i) The ac¬ 
celerated protons and He 2+ may be sharing significant 
energy with each other rather than interacting mainly 
elastically with the background turbulence. If this is the 
case, it will influence all aspects of Fermi acceleration, 
(ii) The particle mean free path may not be a monoton- 
ically increasing function of momentum or it may differ 
substantially for protons and He 2+ . (iii) Different A/Z 
particles may interact differently with the viscous sub¬ 
shock layer. A basic assumption for thermal leakage in¬ 
jection is that the subshock is essentially transparent, 
i.e., phenomena such as cross-shock potentials or large- 
scale turbulence do not strongly influence the injection 
process. If these phenomena are important, it is like¬ 
ly they will influence different A/Z particles differently, 
modifying the A/Z enhancement seen in Fig. [13] 

4. CONCLUSIONS 

As complicated as particle acceleration in relativis¬ 
tic shocks may be, one aspect is profoundly simple: 
if the acceleration is efficient and a sizable fraction 
of the bulk plasma flow energy is put into individu¬ 
al accelera ted particles, as is often assumed in applica¬ 
tions (e.g.. iKulkarni et ~ahllB)fMt iPiran et al.ll2001b iPiranl 
12004 !Meszarosl2006D . the accelerated particles must self- 

consistently modify the shock structure to conserve mo¬ 
mentum and energy regardless of the plasma physics de¬ 
tails. Assuming first-order Fermi acceleration, we have 
investigated how the kinematics of shock modification 
influences the relative acceleration of electrons, protons, 
and heavy elements (i.e., He 2 " 1 ") using a Monte Carlo 
simulation with a dynamic range large enough to model 
acceleration from injection at non-relativistic thermal en¬ 
ergies to ultra-relativistic CR energies. Fig.|4]shows a 12 
decade range in plasma-frame momentum and a greater 
than 20 decade range in dN/dp. A corresponding range 
in photon emission is also obtained (e.g., Fig. hid. 

The underlying wave-particle plasma interactions, 
which are parameterized in the Monte Carlo code, will 
influence details of the shock modification and the resul¬ 
tant radiation; they will determine if acceleration is, in 
fact, efficient and set the maximum energy particles ob¬ 
tain. However, our results show general aspects that are 
largely independent of the poorly known plasma physics 
details if the acceleration is efficient. Considering only 


the kinematics, electrons will be accelerated much less 
efficiently than ions if the shock structure is modified by 
the heavy particles. This result assumes that the heavy 
particles and electrons diffuse in a similar fashion, as 
indicated in equation (JT]) . If this is the case, the A/Z 
enhancement effect we describe increases the injection 
and acceleration efficiency of high A/Z particles com¬ 
pared to low A/Z ones. This dramatically decreases the 
abundance of accelerated electrons compared to heavier 
ions (e.g., Fig. [ 6 ]). The kinematics suggest that relativis¬ 
tic shocks will not be able to place a sizable fraction of 
the shock kinetic energy flux into leptons if protons are 
present. 

Of course, beyond kinematics, the magnetic tur¬ 
bulence produced by wave-particle interactions pla ys a 
critical role and recent re sults ('e.g.. iSironi fe Spitkovskvl 
120111 : iKumar et al.ll2015[ ) show that some fraction of the 
proton energy can be transferred to electrons in the shock 
precursor via magnetic turbulence. These PIC results 
are particularly important for astrophysical applications 
where the radiating electrons presumably contain a siz¬ 
able fraction of the available energy budget. In relativis¬ 
tic shocks, heavy elements must transfer a sizable frac¬ 
tion of their energy to electrons for Fermi acceleration to 
be relevant for electrons. 

We have modeled this energy transfer by including 
a parameter, /; on , that sets the fraction of ion energy 
transferred to electrons as the particles first cross the 
subshock. While the effect of f- Ion is large (e.g., Fig. 0 , 
kinematics must still play a role: light and heavy parti¬ 
cles will be treated differently in relativistic flows. This 
is seen clearly in Figs. [ 6 ] and [11] where, for a given /i on , 
the e/p ratio drops substantially between the UM case 
(where no A/Z effect occurs) and the NL case. Figs. [ 6 ] 
and El where /i on = 0.1 for both the 70 = 10 and 1.5 
shocks, also show that a larger fraction of ion energy 
must be transferred to electrons for high Lorentz factor 
shocks to produce a significant e/p ratio. 

We make a clear prediction that is directly testable 
with PIC simulations. If Fermi acceleration is efficient 
enough so the shock structure is modified by the back¬ 
pressure of accelerated particles, heavy element ions will 
show a clear enhancement over protons (i.e., Fig. H3D- 
We know of no non-kinematic effects (e.g., cross-shock 
potentials, energy transfer via wave-particle interactions, 
or other electrostatic processes) that can produce such 
an enhancement. 

The combined processes of energy transfer from 
heavy particles to electrons, and the kinematics of shock 
smoothing, produce strong signatures on the radiation 
emitted by these particles. In Fig. [T2] we show results 
for 70 = 10 and 70 = 1.5. Since there are a number of 
important parameters that influence the emission, such 
as shock Lorentz factor, ambient density, magnetic field, 
and size of the emitting region, it is non-trivial to char¬ 
acterize the emission. Nevertheless some general prop¬ 
erties stem mainly from the kinematics and should hold 
regardless of the plasma physics details. 

Particle spectra should harden as the shock speed de¬ 
creases from fully relativistic to non-relativistic speeds, 
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mainly because the shock compression ratio increases 
and, for low enou gh 70 , Rtot > -Rrh (see figure 10 in 
lEllison et al.ll201~3h . However, even though the compres¬ 
sion ratio (defined as Rtot = U 0 /U 2 ) is lower for ultra- 
relativistic shocks, the downstream local plasma number 
density n 2 = 7 o/ 3 o«o/( 72 / 32 ) can be large, enhancing the 
pion-decay emission more than IC and synchrotron. The 
possibility of a significant change in the character of non¬ 
linear effects in the trans-relativistic regime, as well as 
the fa ct that trans-relativisti c shocks have been observed 
(e.g., ISoderberg et alJfeOlOf) . makes this an important 
area for future work. 

The magnetic field is a critical parameter for syn¬ 
chrotron emission and, with the exception of Model 
C, we have assumed B 0 = 100 /zG for the background 
field. In our plane-parallel approximation, the back¬ 
ground field remains constant throughout the shock. 
Values of 100’s of /zG can be expected for a shock mov¬ 
ing through an ambient field of a few pG when com¬ 
pression and amplification are considered. Compression 
will increase the field by a factor ~ 70 and NL amplifi¬ 
cation, as believed to occur in strong, non -relativistic 
shocks in young SNRs (see, for examp le, iRdl I 2004 I 
I Vladimirov et a l. 2009: lBvkov et ahll20l4 and references 
therein), may increase the strength further. For simplic¬ 
ity, we have not attempted to included compression or 
magnetic field amplification of t he magnetic fi eld here. 
Field compression is included in IWarrenl (l2015h . 

One important aspect of the changing afterglow 
emission as the shock slows from ultra-relativistic to non- 
relativistic speeds is the position of the synchrotron peak. 
As seen in Fig. [ 12 ] the peak shifts from ~MeV to ^keV 
as the shock slows from 70 = 10 to 70 = 1.5. A detailed 
evolutionary model of GRB afterglows using Monte Car- 
lo techniques for NL Fermi acceleration is presented in 
IWarrenl ( 2015 ). In this afterglow model, the Monte Car¬ 
lo simulation is combined with an analytic or numerical 
description of the jet-shock evolution. The shock accel¬ 
erated particles and resultant radiation are calculated at 
various times as the shock moves through the jet and the 
total emission observed at Earth is determined. 

We caution that, for simplicity, we have assumed 
Bohrn diffusion here (i.e., equation [2] A m f p = r g oc 


p ) whereas the actual scattering process in relativis¬ 
tic plasmas is certain to be more complicated (e.g., 
iLemoine fc Pelletier! 120101) . In particular, particles in¬ 
teracting with the small-scale turbulence generated by 
the Weibel instability are more likely to have A m f p oc p 2 
and additional instabilities may contribute longer scale 
turbulence with a different momentum dependence (e.g., 
iCasse et~aT1l2013l : ILemoine et al .1120 Mil . In fact, the mo¬ 
mentum dependence and normalization of A m f p can be 
expected to vary with momentum as well as position 
relative to the subshock. While a first-principles deter¬ 
mination of A m f p will undoubtedly require PIC simula¬ 
tions, the nonlinear effects we describe here stem from 
basic considerations of momentum and energy conserva¬ 
tion and should persist if first-order Fermi acceleration is 
efficient and A m f p is an increasing function of p. Bearing 
in mind our simple scattering assumptions, we believe 
this work is the first to include electrons, protons, and 
heavier elements in a nonlinear relativistic shock acceler¬ 
ation model. We predict an A/Z enhancement effect for 
heavy ions in relativistic shock acceleration, and include 
the photon emission consistently with nonlinear parti¬ 
cle acceleration. Nevertheless, modifying the scattering 
prescription is certain to produce important quantitative 
differences, particularly for the photon emission which is 
strongly dependent on the maximum CR energy, and we 
are generalizing the Monte Carlo technique to include 
more realistic forms for A m f p which will be presented in 
future work. 
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TABLE 1 

Model Parameters. 


Model a 

Type b 

70 

/ion 

Bo 

iUpS 

±DwS 


Rtot 

Ng 

e H C 

CHe 

Cel 





[/iG] 

r g0 

r g0 







A 

UM 

10 

0 

100 

10 4 


3.02 

3.02 

2000 




B 

UM 

10 

0.15 

100 

10 4 


3.02 

3.02 

2000 




C 

UM 

10 

0.15 

3 

10 3 

10 3 

3.02 

3.02 

2000 




D 

UM 

10 

0.1 

100 

300 

10 3 

3.02 

3.02 

2000 




E 

NL 

10 

0.1 

100 

300 

10 3 

3.02 

3.02 

2000 

0.60 

0.30 

0.10 

F 

NL 

10 

0.4 

100 

300 

10 3 

3.02 

3.02 

2000 

0.36 

0.20 

0.44 

G 

UM 

1.5 

0.1 

100 

300 

10 3 

3.53 

3.53 

200 




H 

NL 

1.5 

0.1 

100 

300 

10 3 

3.53 

3.9 ±0.4 

200 

0.52 

0.27 

0.08 

I 

UM 

30 


100 

300 

10 4 

3.00 

3.00 

10 4 




J 

NL 

30 


100 

300 

10 4 

3.00 

3.00 

10 4 

0.66 

0.34 



a Models A, B, and C have n p — n e = 1cm 3 with no helium. Models D—H have n p = 1cm 3 , nn e — 0.1cm 3 , and n e = 1.2cm 3 . 
Models I and J have n p = 1cm -3 and nn e = 0.1n p without electrons. The far upstream temperature is 10 6 K in all cases and all models 
have 77 m f p = 1. 

b In the nonlinear (NL) models the shock structure is determined self-consistently. The unmodified (UM) models have a discontinuous 
shock structure with no shock smoothing. 

c For self-consistent NL models, this is the fraction of total energy placed in all particles with energies above 100 MeV as measured in the 
shock frame. 






